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1 .  Introduction 


Multi-dimensional,  time-dependent  flows  are  generally 
considered  as  problems  which  can  be  analyzed  numerically  only  by 
high-speed,  large-scale  computers,  because  of  the  large  number  of 
grid  points  necessary  to  achieve  the  required  degree  of  accuracy. 
The  need  is  emphasized  when  viscous  effects  have  to  be  taken  into 
account  [1].  There  is  no  doubt,  however,  that  the  inclusion  of 
viscous  effects  is  a  must  to  achieve  physical  reality  in  certain 
problems.  To  this  category  belong  all  flows  involving  separa¬ 
tion,  plume  formation,  shock-boundary  layer  interaction  with 
upstream  propagation  of  signals  in  a  generally  supersonic  flow, 
self-sustained  flutter,  etc. 

The  use  of  appropriate  numerical  techniques  reduces  the 
need  for  a  large  number  of  grid  points  and  allows  a  mini-computer 
to  be  used  for  the  analysis.  After  experimenting  with  a  certain 
number  of  such  problems  and  finding  that  the  above  statement  can 
be  supported  by  concrete  evidence,  I  consider  proper  to  report 
the  techniques  in  full  detail.  Two-dimensional  o»-  axi-symmetric , 
time-dependent,  viscous  flow  problems  at  high  Reynolds  numbers 
will  be  considered. 

We  will  discuss  the  basic  points  which  make  the  technique 
efficient,  viz.: 

a)  a  proper  formulation  of  the  equations  of  motion, 

b)  the  use  of  mappings, 

c)  the  stretching  of  coordinates, 

d)  the  discretization  of  the  equations, 

e)  the  treatment  of  imbedded  shocks. 


2.  Equations  of  motion 


In  vector  form,  the  equations  of  motion  for  a  viscous 
flow  (Navier-Stokes  equations),  assuming  that  the  viscosity,  u  , 
i3  a  constant,  are: 
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Equations  of  motion 
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where  p,  p,  S  and  0  are  the  thermodynamical  parameters  density, 

pressure,  entropy  and  temperature,  respectively.  In  the  same 
♦ 

equations,  V  is  the  velocity  vector,  t  is  the  time,  <  is  the 
coefficient  of  heat  conduction,  and  *  is  the  dissipation  term. 
The  latter  is  a  well-known  non-negative  quadratic  form  depending 
on  the  space  derivatives  of  the  velocity  components.  Its  expres¬ 
sions  for  the  various  orthogonal  coordinate  systems  considered 
in  the  present  Report  will  be  given  later  on.  If  pressures,  den¬ 
sities  and  lengths  are  expressed  in  terms  of  reference  values, 

p  ,  p  and  x  ,  respectively,  corresponding  units  of  veloci- 
ref  ref  ref 

ty,  time  and  temperature  are  defined  as 


u  =  (p  Jo  J 

ref  ref  ref 


1/2 


t  =  x  /  u 
ref  ref  ref 


0  ,  =  u  J R  (2) 
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where  R  is  the  gas  constant. 


A  Reynolds  number  and  a  Prandtl  number  can  be  defined  as 


R 
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o  J  Jv, 
ref  ref  ref 


P  =  c  \i/k 
r  p 


where  c  is  the  specific  heat  at  constant  pressure.  Finally,  let 

the  uni?  entropy  be  c  (the  specific  heat  at  constant  volume).  The 

v 

equations  of  motion  can  thus  be  written  in  the  form: 
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In  an  inviscid  flow,  the  terms  here  affected  by  1/R  do 

e 

not  appear.  For  a  proper  numerical  analysis  of  such  flows,  whose 
mechanics  is  governed  by  the  propagation  of  sound  waves  and  by 
the  convection  of  entropy  along  particle  paths,  it  is  convenient 
to  consider  the  logarithm  of  pressure,  P,  instead  of  the  density, 
as  an  unknown,  and  to  make  use  of  the  equation  of  state  for  a 
perfect  gas: 


3  =  y  In  9  ~  (y-1  )P 


(4) 


so  that  the  equations  of  motion  for  an  inviscid  flow  (Euler's 
equations)  can  be  written  in  the  form: 
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(5) 


If  the  flow  is  viscous,  the  basic  phenomena  of  wave  pro¬ 
pagation  are  still  present,  although  modified  by  the  concurrent 
effects  of  diffusion.  Therefore,  it  would  not  be  advisable  to 
drop  the  basic  integration  techniques  for  convective  terms;  we 
will  consequently  write  the  Navier-Stokes  equations  in  the  form: 
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(where  the  material  derivatives  have  been  replaced  by  partial 
derivatives,  as  costumary,  and  q  is  the  modulus  of  the  velocity). 
In  what  follows,  the  density  will  no  longer  be  used  explicitly, 
and  p  will  be  redefined  in  Eq.  (8). 
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3_.  Conformal  mapping 


We  will  consider  two  types  of  flows,  both  depending  on 
two  space  variables:  two-dimensional  (plane)  flows  and  axi- 
symmetric  flows.  Two  Cartesian  coordinates,  x  and  y,  will  be 
used  for  two-dimensional  flows,  the  same  symbols  will  be  used 
for  the  axial  and  radial  coordinate,  respectively,  in  any  meri¬ 
dional  plane  of  an  axi-symmetric  flow.  The  (x.y)-plane  will  be 
called  the  physical  plane. 

In  addition,  we  will  introduce  a  complex  variable 


z  =  x  +  i  y 


(7) 


and  assume  that,  in  general,  the  physical  plane  will  be  confor¬ 
mally  mapped  onto  an  auxiliary  plane,  described  in  terms  of  a 
complex  variable  c. 


Obviously,  if  no  mapping  is  needed,  all  the  formulas  below  are 
valid,  provided  that  ;  =  z,  that  is,  5  =  x  and  n  =  y.  The  map¬ 
ping  function  must  be  chosen  in  such  a  way  that  the  contours  of 
the  flow  field  are  as  close  as  possible  to  £=constant  and 
n=constant  lines,  or  to  p=constant  and  9=constant  lines  [2].  Let 


d?  iui 
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Let 


We  will  begin  with  the  case  in  which  £  and  n  are  used. 


1_  d  log  g 
g  dz 


+  i  ♦, 


(12) 


Let  also  I ,  3  be  the  unit  vectors  tangent  to  the  n  =  constant 
line  and  to  the  5  =  constant  line  in  the  z-plane,  and  u  and  v  be 
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Conformal  mapping 

N 

the  corresponding  velocity  components,  respectively. 

so  that 

t 

V  =  u  I  +  v  3 

(13) 

It  is  convenient  to  note  that 

£  =GZ ,  £  =GS  •  n  s— GS ,  ri  •GZ 
x  y  x  y 

(14) 

x  =Z/G,  X  =-3/G,  y  =3/G,  y  =0/G 

5  n  5  n 

(15) 

2  2  2  2  2  2 

5  +5  =G  ,  x  +x  =1/G 

x  y  K  h 

(16) 

/ 

G  =G$ ,  G  =— G<t>  ,  id  =*  ,  u>  =$ 

C  1  n  2  5  2  n  1 

(17) 

If  I  and  J  are  the  unit  vectors  parallel  to  the  x- 
respectively ,  and 

and  y-axis. 

V  =  U  I  ♦  V  J 

> 

(18) 

it  follows  that 

i.i=e  ,  £.3*8  .  J.I=-3  ,  3-3=e 

(19) 

U  =u^-v2,  V=u 3+v£,  uzuZ+V3,  v=-U3+V0 

(20) 

For  any  element  ds  in  space,  we  have 

2  2  2  2  2  2 
ds  =  d£  /G  +dn  /G  +dx^ 

(21) 

if  x  is  the  third  Cartesian  coordinate,  in  a  two 
problem,  and 

-dimensional 

2  2  2  2  2  2  2 
ds  =  d£  /G  +dn  /G  +y  dx 

3 

(22) 

r 

if  x  is  the  angular  coordinate  in  an  axi-symmetric  problem. 

Consequently,  in  a  two-dimensional  problem,  the 

basic  vector 

operators  in  ( 1 )  can  be  expressed  as  follows: 

VP=G(P  1+P  j) 

5  n 

(23) 
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•»  2  v  u 
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(24) 

VxVxV=G2[(^)  -£)  ](vi-uj) 

G  5  G  n 

(25) 
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(26) 
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55  nn 

(27) 

2  2  2  2  2  2  2 
$  =  4(e  +e  +e  )  +  "^[(e  -e  )  +(e  -e  )  +  (e  -e  )  ]  (28) 
12  23  31  3  11  22  22  33  33  11 
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e,  =G(u  +v<t  ),  e  =G(v  -u<t>  ),  e  =  —  [(Gv)  +(Gu)  ]  (29) 

11  52  22  nl  12  2  5  n 

and  e^,  e^,  e^  equal  to  zero. 

In  an  axi-syrametric  problem,  (23),  (24),  (25)  and  (28) 

still  hold,  but  (26)  and  (27)  must  be  replaced  by 

*  2  u  v  V 

7 . V  =  G  [(-)  +(“)  ]  +  -  (30) 

G  5  G  n  y 

2  2  G 

7  9  =G  (9  +9  )  +  (39  +29  )  (31  ) 

55  nn  V  5  n 

and  e ^  is  not  zero,  but 

V 

e  =  ~  (32) 

33  y 

From  now  on,  we  will  follow  the  common  practice  of  using  a  multi¬ 
plier,  j,  which  is  equal  to  zero  for  two-dimensional  flows  and  to 
1  for  axi-symmetric  flows.  We  will  denote  by  a  the  only  non-zero 
component  of  the  curl  of  V  (24): 


fl  =  G  (v  -u  —  v<4  — u*  ) 
5  n  1  2 

and  by  A  the  divergence  of  V  (26)  or  (30): 


1  =  *1.^22^33-  *33  *  J  7 

and  introduce  the  symbols 
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D  =  G  ( v$  i  +u$  ) 


E  s  G(v$  — )  +  e 

2  1  33 


(35) 


It  is  easy  to  see  that,  in  terms  of  the  independent  variables  £ 
and  n,  the  Navier-Stokes  equations  are: 


P  +  G(uP  +  vP  )  +  yG(u  +  v  )  +  yE  =  — 
t  5  n  ?  n  Dt 

4  l  ,  9 

u  +  G(uu  +  vu  )  +  vD  +  G9P  =  [-GA  -Gfi  -j  —  n]  — — 
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4  2  e 

G(uv  +  vv  )  -  uD  +  GeP  =  [—  Ga  +Gq  +j  —  n]  - 
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S  +  G(uS  +  vS  )  =  C  (y-1  720  ]/( pR  ) 

t  5  n  P  e 
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(36) 


with  $  defined  by  (28)  and 

2  2  G 

7  9  =  G  (9  +9  )+j  -  (29  +Z9  )  (37) 

((  y  C  n 

If  the  basic  variables  in  the  c-plane  are  p  and  9 ,  we 
will  define  i  and  j  as  the  unit  vectors  tangent  to  the  9=constant 
lines  and  p=constant  lines,  respectively,  and  the  velocity  com¬ 
ponents,  u  and  v,  will  be  defined  accordingly,  so  that  (13)  still 
holds.  Practically,  all  equations  from  (12)  through  (37)  must  be 
changed,  as  follows. 


_  £  d  log  g 

♦  '  g  dz  '  *1  +  42 

(38) 

G  c  1  ( 9  -o> ) 

Z+i2  =  —  =  e 

Pg 

(39) 
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(40) 

G 
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=  pZ/G 

(42) 

2  2  2  22,2  ,2 
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(43) 
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v  G 
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1  V  1 
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V 
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33  y 

(50) 

G 

D  =  — C  —  v(  1  — ,  )+u$  ] 

P  1  2 

(51) 

G 

E  =  — Cu(  1  — 4> ,  )-*-v<t>  ]  +  e 

P  1  2  33 

(52) 

1 

Q  =  G(v  -~u  )  -  D 
p  p  9 

(53) 

A  =  e  +  e  +  e 

11  22  33 

(54) 

In  this  case,  the  Navier-Stokes  equations  become: 


v  8  DS 

P  +  G(uP  ♦  P  )  +  fG(u  +  )  +  y E  = 

t  p  p  8  p  p  Dt 

v  4  G  e  9 

u  G(uu  +—  u)+vD+G0P  =  (~Ga  -  —  fl  -  j  TJ)  - 

t  p  P  9  p  3  p  P  9  y  pR 
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v  G  4  G  Z  Q 

v  +  G(uv  +~v)-uD+~9P  =  (t  ~  A  +Gq  +  j  ~n)  - 

t  ppe  p  9  3  p  9  p  ypR 
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S  +  G(uS  +  -S  )  =  [(Y-1)*  +  ^-17  9]/(pR  ) 
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(55) 
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with  *  defined  by  (28)  and 

2  9  2+9  2/p 

7^9  =  —  (9  +p9  +~  9  )  +  jG  — - -  (56) 

p  p  pp  p  99  y 

4^  Computational  plane 

The  c -plane,  however,  is  not  the  computational  plane.  In 
most  instances  (and  particularly  for  viscous  flows,  when  boundary 
layers  and  shear  layers  are  present),  grid  lines  must  concentrate 
in  the  regions  of  highest  gradients.  An  efficient  way  for 
achieving  an  uneven  partition  of  grid  lines  consists  of  using 
stretching  functions,  which  essentially  define  a  correspondence 
between  n  (or  p)  and  a  new  variable,  X,  and  between  5  (or  9)  and 
a  new  variable,  Y.  The  boundaries  of  the  region  to  be  computed 
in  the  physical  plane  should  be  made  to  correspond  to  the  X=0, 
X=1  lines  and  to  the  Y=0,  Y=1  lines,  respectively.  The  computa¬ 
tion  is  performed  in  the  (X,  Y)  plane,  using  a  rectangular  grid 
with  evenly  spaced  lines;  it  must  be  noted,  however,  that 
X=constant  lines  and  Y=constant  lines  on  the  physical  plane  are 
generally  not  orthogonal  to  each  other,  since  the  boundaries  are 
not  necessarily  represented  by  C=constant  lines  and  n=constant 
lines,  or  by  psconstant  lines  and  9=constant  lines. 

Wc  will  consider  here  two  examples.  In  the  first,  the 
basic  variables  in  the  ;-plane  are  £  and  n;  two  boundaries  are 
defined  by  constant  values  of  ?  and  two  boundaries  are  defined  by 
values  of  n  which  depend  on  £  and  t,  as  well.  It  will  be  neces¬ 
sary,  thus,  to  define  Y  not  only  as  a  function  of  n  but  of  £  and 
t,  and  to  distinguish  carefully  between  the  time,  t,  in  the  phy¬ 
sical  plane  (or  in  the  t-plane)  and  the  time,  T,  in  the  computa¬ 
tional  plane.  Therefore,  we  will  write: 

X  =  X(5) 

Y  =  YU.n.t)  (57) 

T  =  t 
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Letting 


a  =  GuX  ,  a  =  yGX  ,  b  =  Y  +GuY  +GvY  ,  b  =  yGY 

11  5  12  5  11  t  5  n  12  ^8) 

g  s  yGY  ,  a  -  GeX  ,  h  s  GsY  ,  b  =  G0Y 

11  5  21  c  2  5  21  n 


F  =  [  ( y —  1  )*  +  729]/(pR  ) 

P  e 


ci  =  yE  -  F 

g  4  £ 

o  =  vD  -  ~ —  [-GA  -Gfl  -j  _  £2] 

2  pR  3  5  n  y  .... 

e  (60) 

0  4  2 

c  =  -  uD  -  - [“Ga  +G£2  +j  -  a] 

3  pR  3  n  Z  y 

e 

the  equations  of  motion  are  recast  in  the  form: 

P  +a  P  +b  P  +a  u  +g  u  +b  v  +  c  =0 

T  11  X  11  Y  12  X  11  Y  12  Y  1 

u  +a  u  +b  u  +a  P  +h  P  +c  =0 
T  11  X  11  Y  21  X  2  Y  2 

(61 ) 

v  +a  v  +b  v  +b  P  +c  =0 
T  11  X  11  Y  21  Y  3 

VhWy =  F 

In  the  second  example,  the  basic  variables  in  the  c-plane 
are  p  and  e;  two  boundaries  are  defined  by  constant  values  of  9 

and  two  boundaries  are  defined  by  values  of  p  which  depend  on  9 

and  t.  It  will  be  necessary,  thus,  to  define  X  not  only  as  a 

function  of  p  but  of  9  and  t;  therefore,  we  will  write: 

X  =  X(p ,9 ,t) 

Y  =  Y (9 )  (62) 
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Letting 


a  =  X  +GuX  +GvX  /p  ,  a  =  yGX  ,  b  =  GvY  /p 
11  t  p  0  12  p  11  0 


b12  =  yGY0/p  1  g12  =  YGX0/P  *  a21  =  G9Xp 


(63) 


h  =  GeX  /p  ,  b  =  GeY  /p 
3  0  21  0 


and  F,  c  as  above  in  (59)  and  (60),  respectively,  and 

0GJ  ,  Z  , 

c  =  vD  -  -—[-a  -n  /p  -  j  r-n] 

2  pR  3  p  0  Gy 

e 

„  0G  4  .  8  _ 

C  =  -uD  -  - [-A  / p+Q  +  J 

3  pR  3  9  P  Gy 

e 


(64) 


the  equations  of  motion  are  recast  in  the  form: 


P 

T 


+a  P  +b  P  +a  u  +b  v  +g  v  +c 
11  X  11  Y  12  X  12  Y  B1 2  X  1 

u  +a  u  +b  u  +a  P  +c  =0 
T  11  X  11  Y  21  X  2 


v  +a  v  +b  v  +b  P  +h  P  +c  =0 
T  11  X  11  Y  21  Y  3  X  3 


S  +a  S  +b  S  =  F 
T  11  X  11  Y 


=0 


(65) 


It  is  easy  to  see  that  the  first  three  equations  (61)  and  the 
first  three  equations  (65)  have  the  same  form  as  Equations  (23) 
of  [3].  The  integration  procedure  explained  in  [3]  can  be  ap¬ 
plied.  We  expect  the  integration  technique  to  provide  a  very 
good  estimate  of  convection  and  wave-propagation  effects  which 
are  common  to  inviscid  and  viscous  flows. 
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5..  Rigid  boundaries 


For  a  viscous  flow,  the  velocity  at  a  rigid  boundary  is 
assumed  to  vanish.  There  are  no  difficulties,  thus,  in  the 
determination  of  u  and  v.  Pressure  and  entropy  require  more 
care.  We  consider  here  the  case  of  an  isothermal  wall;  the  tem¬ 
perature,  9  is  a  prescribed  constant.  Because  of  (4),  P  and  S  are 
linearly  related;  therefore. 


DS 

DT 


O-y) 


DP 

DT 


and  consequently,  the  first  of  (36)  and  the  first  of  (55)  must  be 
replaced  by 


and 


P  +  G(uP  +vP  )  +  G(u  +v  )  =  0  (66) 

*  5  n  5  n 

P  +  G  (uP  +~P)+G(u  +  -  v  )  =  0  (67) 

t  p  p  9  p  p  9 


respectively,  having  taken  into  account  that,  in  both  cases,  E 
vanishes  identically.  The  equations  to  be  solved  are  still  (61.) 
and  (65),  provided  that  a^,  b^,  and  g^  are  divided  by  y 

and  c^  is  set  equal  to  zero. 


Once  the  pressure  has  been  determined,  S  is  obtained  from 
(4).  Therefore,  S  is  not  calculated  directly  as  a  result  of  dis¬ 
sipation  and  heat  transfer  in  the  flow,  because  the  condition  of 
constant  wall  temperature  implies  some  external  action,  and  (4) 
gives  us  the  final  outcome  of  such  action  and  the  local  variation 
of  pressure. 


S.  Integration  procedure 


The  equations  of  motion  (61)  or  (65)  are  integrated  fol¬ 
lowing  the  general  guidelines  of  Section  6  in  [3].  We  define 


X 

C  sg  u  +  c 
1  1 1  Y  1 


X 

C  =b  u,  ♦  h  P  c 
2  11  Y  2  Y  2 


12 


I 


Integration  procedure 


„Y 

C  =g  v 
1  12  X 


,  C  :a  v  +  h  P  +  c 
2  11  X  3  X  3 


X  Y 

After  finding  the  characteristic  slopes,  X  ,  X  (i=1,2)  from 

l  i 


X 

a  -X 

a 

Y 

b  -X 

b 

11  i 

21 

11  i 

21 

a 

X 

a  -X 

=0  , 

b 

Y 

b  -X 

12 

11  i 

12 

11  i 

and  letting 


X  Y 

a  -X  b  -X 

11  i  11  i 

A  =  -  ,  B  =  - 

i  X  X  i  Y  Y 
X  -X  X  -X 

2  1  2  1 


a  b 

DX  a—ii-  dY  - 

ij  XX’  ij  Y  Y 

X  -X  X  -X 

2  1  2  1 


the  equations  to  be  integrated  are: 

X  .  X„  X„  X  ,  X  X  .  X  „ 

PT  +  A1X1PX1_A2X2PX2  +  °12  X2UX2_X1UX1  +  C1 

UT  +  °21 (X2PX2_X1PX1 )  +  A1X2UX2_A2X1UX1  +  °2  =° 


Y  Y  Y  Y  Y  Y  Y 

P  +BXP  -BXP  +  D  (Xv  -Xv  )+C  =0 

T  1  1  Y1  2  2  Y2  12  2  Y2  1  Y1  1 

Y  Y  Y  Y  Y  Y 

v  +D  (XP  -XP  )  >  B  X  v  -B  X  v  +C  =0 
T  21  2  Y2  1  Y1  1  2  Y2  2  1  Y1  2 


X  Y 

PT  =  PT  ♦  PT  (73) 

TIT 

For  discretization,  the  space  derivatives  are  classified 
into  three  categories:  (i)  the  ones  explicitly  appearing  in  (71) 
and  (72),  (ii)  the  ones  explicitly  appearing  in  (68),  and  (iii) 
the  ones  appearing  in  (59),  (60)  and  (64).  The  derivatives  of 
the  first  category  are  discretized  according  to  the  rules  (14) 
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and  (15)  of  [31;  the  derivatives  of  the  second  category  according 
to  the  rules  (34)  and  (35)  of  [33;  an  the  derivatives  of  the 
third  category  are  approximated  by  ordinary  centered  differences, 
because  the  physical  nature  of  the  terms  which  they  affect  is 
diffusion.  Few  exceptions  to  the  general  rules  are  necessary  for 
points  on  rigid  boundaries  or  next  to  rigid  boundaries.  For 
points  on  rigid  boundaries,  the  alternate  two-point-three-point 
approximation  is  always  taken  using  points  inside  the  flow  field. 
For  points  next  to  rigid  boundaries,  if  use  of  a  point  located 
behind  the  wall  were  required  in  a  three-point  approximation,  the 
latter  is  substituted  by  a  two-point  formula. 


7.  Numerical  treatment  of  shocks 


Shocks  are  generally  present  in  a  compressible  flow 
field,  either  as  boundaries  or  imbedded  in  the  flow.  In  both 
cases,  they  can  be  treated  numerically  in  the  general  framework 
of  the  computational  technique  described  in  the  previous  Sec¬ 
tions.  We  begin  with  some  general  considerations. 

Let  N  and  x  be  the  unit  vector  normal  and  tangential  to  a 
shock  at  any  of  its  points,  Q,  respectively: 

*  *  V  *  M23  ’  ?  =  “V  +  M13  (74) 

W  the  shock  velocity, 

W  =  W  R  (75) 

and  u,  v  the  velocity  components  along  fl  and  x,  respectively: 

u  =  uN  +vN  u  =  uN  -vN 

12  _  1  „  2  (76) 

v  s  -uN  +vN  v  =  uN  -t-vN 

2  1  2  1 

The  N-component  of  the  velocity  relative  to  the  shock  is 
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u  =  U  -  W  (77) 

rel 


The  relative  normal  Mach  number  on  the  low  pressure  side  of  the 
shock  is 


„2 

nlrel  yei 


(78) 


The  Rankine-Hugoniot  conditions  are: 


2yM  t  ,  -Y+1 
nlrel 

P  =  P  +  In - 

2  1  Y+1 


y9 

f-1  ,  21 

u .  + 


(79) 


2rel  y+1  Irel  y+1  u 


Irel 


to  W: 


We  will  need  the  derivatives  of  P  and  u  with  respect 

2  2rel 


3  P 


3W 


4u 


Irel 


<er(*-1)9i 


(80) 


3U 


2rel 


3W 


Ye 


Y  +  1  ~2 

U 

Irel 


Let  us  assume  that  the  shock  is  oriented  in  the  general 

direction  of  the  X=constant  lines;  therefore,  it  can  be  defined 

by  its  intersections,  X  ,  with  Y=constant  lines.  At  each  inter- 

s 

section,  we  consider  two  points,  one  on  the  low-pressure  side  and 
the  other  on  the  high-pressure  side.  The  point  on  the  low- 
pressure  side  must  be  updated  by  using  information  proceeding 
from  that  side  only.  The  point  on  the  high-pressure  side  must  be 
updated  by  using  information  from  both  sides.  The  information 
proceeding  from  the  low-pressure  side  must  satisfy  the  Rankine- 
Hugoniot  conditions;  the  information  proceeding  from  the  high- 
pressure  side  is  carried  to  the  shock  along  a  characteristic. 
The  acceleration  of  the  shock  results  from  the  compatibility  of 
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the  different  information.  To  obtain  the  characteristic  equa¬ 
tion,  we  should  rewrite  the  equations  of  motion  in  a  frame  rela¬ 
tive  to  the  moving  3hock,  where  the  Y-coordinate  of  the  shock  is 
unchanged  but  its  X-coordinate  follows  the  shock  in  its  motion. 
Therefore,  we  will  introduce  a  new  set  of  coordinates, 


X  =  X  -  X  (Y.T) 
s 


e  =  Y 


(81) 


x  =  T 


and  rewrite  (61)  and  (65)  in  the  form: 


P  +a  P  +b  P  +a  u  +b  v  +y  v  +c  =0 
t  11  x  He  12  x  12  e  12  X  1 

u  +a  u  +b  u  +a  P  +h  P  +c  =0 
t  11  x  lie  21  x  2e  2 

v  +<tM  v  +b  v  +b  P  +n  P  +c  =0 
t  11  X  11s  21s  3  x  3 


(82) 


S  +a  S  +b  S  =F 
x  11  x  He 


where 


ai1=ail"b11XsY~XsT  ’  ai2=ai2-g1 1XsY  ’  a21=a2rh2XsY 
Y12“g12~b12XsY  ’  n3=h3-b21XsY 


(83) 


A  characteristic  equation  may  now  be  obtained,  using  x  and  t  as 
basic  independent  variables: 


(a  -X)(P  +XP  )— a  (u  +Xu  ) — y . _ ( v  +Xv  )+R  =  0  (8^) 

11  t  x  12  t  x  12  x  x 

where  X  is  the  solution  of  the  equation 


a11~x 

“21 

n3 

°12 

a  -X 

11 

0 

=  0 

(85) 

o 

(M 

>- 

ail’X 
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that  is. 


X=a  ±(a  a  +Y  n  ) 
11  12  21  12  3 

and  R  comprises  all  the  remaining  terms. 


Note  now  that 


a  =  yG''N1  .  y12  =  yGvN2 

°21  =  G0vNi  *  n3  =  G9vN2 


where 


2  2  1/2 
(ai2+T12) 


Therefore , 


ai2U  +  Y12V  =  yGvS 


a  u  +  v  v  =  yGvu  -yGv(uN  +vN  ) 
12  x  12  x  x  lx  2x 

a  u  +  y  v  =  yGvu 

12  X  12  x  X 


Instead  of  ( 84 )  we  can  write: 


(o  -X)(P  +XP  )-yGv(u  +Xu  )+R,  =  0 

11  x  x  x  x  1 


where 


R  =  R+yGv (uN  -vN  ) 

1  lx  2x 

A  further  simplification  is  obtained  by  observing  that 

1/2 

.  -H-.±  12  21.— 3i  i -  !ti  (91) 

tGv  yGv  y 

because  of  (87)  and  (88),  so  that  (90)  can  be  written  in  the 
form: 
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±r-  (P  +XP  )+u  +Xu  +R  =  0  (92) 

Y  t  x  t  x  1 

The  values  of  P  and  u  ,  obtained  by  integrating  the 

T  T 

Navier-Stokes  equations  (82)  using  information  from  the  high- 
pressure  side  of  the  shock  only,  must  satisfy  (92).  On  the  other 
hand,  the  exact  solution  of  the  flow  problem  in  the  presence  of 
the  shock,  which  accounts  also  for  the  information  from  the  low- 
pressure  side  and  the  Rankine-Hugoniot  conditions,  must  satisfy 
(92)  as  well.  Therefore,  if  we  write  (92)  twice,  once  for  the 
t-der ivatives  as  obtained  from  the  Navier-Stokes  equations  (that 
is,  using  for  the  shock  point  the  same  integration  procedure 
which  is  applied  to  other  points)  and  again  for  the  exact 
t-der ivatives  and  we  subtract  one  equation  from  the  other,  the 
simple  result  is  obtained: 


a  (NS)  .  .(NS)  n 

:  (P  -P  )  +  u  -u  =  0 

Y  T  T  T  T 


(93) 


where  the  derivatives  obtained  from  the  Navier-Stokes  equations 
are  labeled  (NS)  and  the  exact  derivatives  are  unlabeled. 


The  latter  derivatives  can  obviously  be  expressed  in  the 


form : 


3P2 

P  =  P»  +  -  W 

T  T  3W  T 


3U 2 

u  =  u*  +  -  W 

T  T  3W  T 


(9«) 


where  f*  is  a  derivative  computed  considering  W  as  a  constant. 
The  acceleration  of  the  shock  is  thus  obtained: 


/n(NS)  D_.  ..(NS) 

±a(P  -P»)+y(u  -u*) 

^  _ T _ T _ T _ T _ 

t  ±a  3  P  /3W  +  y  3 u  /3W 


(95) 


Since  both  the  starred  derivatives  and  the  derivatives  indicated 
by  (NS)  are  computed  using  the  same  initial  values,  (95)  can  be 
replaced  by 


W  4t 


(NS )  ..(NS ) 

±a(P  -P*)-*nr(u  -u*) 

±a  3?  /3 W  y  3u./3W 
2  2 


(96) 
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where  P*  and  u*  are  the  values  on  the  high-pressure  side  of  the 
shock  obtained  from  updated  values  on  the  low-pressure  side  by 
applying  the  Rankine-Hugoniot  conditions  to  a  shock  whose 
geometry  has  been  updated  but  whose  velocities,  W,  are  still  the 
same  as  at  the  beginning  of  the  integration  step.  This  result 
[4]  is  remarkable  since  it  provides  an  extremely  simple  method 
for  calculating  the  shock  acceleration  although  it  relies  on  the 
3ame  basic  concepts  which  have  been  shown  to  be  necessary  for  a 
proper,  physically  correct,  handling  of  shocks  [5]. 


We  will  now  consider  first  the  case  where  the  shock  is  a 
boundary,  mcving  into  a  gas  at  rest,  and  then  the  case  of  an  im¬ 
bedded  shock. 


In  the  case  of  a  shock  moving  into  a  gas  at  rest,  let  us 
assume  that  the  3hock  is  a  right  boundary  of  the  computational 
field.  It  is,  thus,  defined  by  X=1 ;  in  this  case,  x.  e  and  x 
coincide  with  X,  Y  and 


4 12_a  12*  Y12_g12’  n3’h3’ 
from  the  high-pressure 


T,  respectively  and  a  =a  ,  a-.|=a  , 

The  characteristic  reaching  the  shock 


side  is  a  right-running  characteristic, 
and  in  the  preceding  equations,  whereas  a  ±  appears,  the  upper 
sign  must 
fore,  u. 


be  used. ^  Note  ^n  addition  that  u^=0  and  v^sO;  there- 

=-W  and  M  =W  /y.  In  this  case,  obviously,  the 
Irel  nlrel 

low-pressure  side  values  are  known  without  any  need  for  comput¬ 
ing;  the  (NS)-values  on  the  high-pressure  side  are  obtained  to¬ 
gether  with  and  using  the  same  procedure  as  for  interior  grid 
points. 


For  an  imbedded  shock,  whose  location  does  not  generally 
coincide  with  a  grid  line,  we  use  a  simplified  procedure  to  ob¬ 
tain  the  values  at  the  shock  on  the  low-pressure  side  and  the 
(NS)  values  on  the  high-pressure  side.  On  the  low-pressure  side, 
instead  of  integrating  (82)  (which  would  require  a  special,  and 
not  easy,  redefinition  of  approximations  for  the  x-  and 
e-derivatives),  we  simply  assume  that  the  values  at  the  shock  can 
be  extrapolated  from  ohe  two  adjacent  grid  points  on  the  same 
ysconstant  line,  both  at  the  end  of  the  predictor  and  the  correc¬ 
tor  level.  On  the  high-pre3sure  side,  we  assume  that  the  T- 
derivatives  on  the  shock  are  equal  to  the  T-derivatives  at  the 
grid  point  next  to  the  shock  on  the  same  Y=constant  line.  Note, 
however ,  that 
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f 

T 


f  +  f  X  „ 
T  x  sT 


(97) 


for  any  function  f.  The  values  of  f  at  the  shock  are  assumed  to 

T 

be  the  same  as  at  the  next  grid  point  on  the  hight-pressure  side, 
on  the  same  Y=constant  line.  The  values  of  f  are  approximated 

X 

as  follows: 


f  * 
X 


[f  -f 
A  s 


+  - -  (f  -f  )]/AX 

1+e  B  s 


(98) 


where  e=(X  -X  )/AX  and  A  is  the  grid  point  next  to  the  shock  on 
A  3 

the  high-pressure  side  and  B  is  the  next  grid  point.  This  formu¬ 
la  provides  a  smooth  transition  when  the  shock  crosses  an 
X=constant  line. 


For  an  imbedded  shock,  thus,  the  calculation  proceeds  as 
follows:  In  the  predictor  stage,  after  updating  all  grid  points, 
the  low-pressure  side  of  the  shock  is  obtained  by  extrapolation 
and  (97)  is  also  applied  to  P,  u,  v  and  S.  The  values  on  the 
high-pressure  side  are  updated  by  adding  f  At  to  th^^nitial 
values  of  f;  the  values  so  obtained  are  the  predicted  f  .  The 
predicted  f*  are  obtained  by  applying  (79)  to  the  predicted 
values  on  the  low-pressure  side.  Then,  (96)  is  applied  and  W  is 
temporarily  updated,  but  its  original  value  is  retained  in 
storage.  The  geometry  of  the  shock  is  updated,  considering  that, 
in  virtue  of  (75),  (14)  and  (19), 

C  .  =  G  W/N  (99) 

st  1 

or,  using  (41)  in  lieu  of  (14), 

o  =  G  W/N  (100) 

St  1 

and  using  the  approximations: 

1  2 

5  (t+At)  =  5  (t)+€  At+-$  At 

s  s  st  2  stt 

1  2  (101) 

p  (t+At)  =  p  (t)+p  At+~ p  At 

3  3  St  2  Stt 

where  the  second  derivatives  are  obtained  by  differentiating  (99) 
or  (100).  In  the  corrector  stage,  the  procedure  outlined  above 
for  the  predictor  stage  is  repeated  through  the  application  of 
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(97).  The  updating  of  the  values  on  the  high-pressure  side  is 
obtained  by  adding 


f(pred) 

x 


)  At 


to  the  predicted  values.  The  corrected  f*  are  obtained  by  apply¬ 
ing  (79)  to  the  corrected  values  on  the  low-pressure  side,  with 
the  original  value  of  W.  Then,  (96)  is  applied  again  and  W  is 
definitively  updated.  The  Rankine-Hugoniot  conditions  (79)  are 
applied  once  more  to  the  corrected  values  on  the  low-pressure 
side  using  the  final  value  of  W,  to  obtain  the  final  values  on 
the  high-pressure  side.  The  entropy  is  also  computed  from 


S  =  S  +  P-P-  yln(u  /u  ) 
2  1  21  Irel  2rel 


(102) 


At  this  3tage,  it  is  convenient  to  correct  the  values  at 
the  grid  point  next  to  the  3hock  on  the  high-pressure  side. 
Pressure  and  velocity  components  are  interpolated  from  the  values 
at  the  shock  and  the  following  grid  point.  Entropy  is  also  in¬ 
terpolated  considering  that  it  is  carried  along  a  streamline. 


8.  An  example 


To  illustrate  the  technique,  let  us  consider  the  follow¬ 
ing  example.  A  duct  (either  two-dimensional  or  axisymmetric)  has 

its  centerline  on  the  x-axis;  its  lower  wall  (AC)  is  defined  as 

the  image,  in  the  z-plane,  of  a  straight  line,  8=9  ,  in  the 

o 

;-plane  (Fig.  1).  The  mapping  is  defined  by 

z  =  (r  /t)  [  (z^-1/z^)/2  -  logz^  -  iir] 

0  2  2  2  (103) 

z2  +  1/z2  =  2B  u  +  1/c) 

where  r  is  an  arbitrary  parameter,  and  B  is  determined  to  assure 
o 

correspondence  between  c=1  and  z=-i.  The  section  of  the  duct  at 
x=0  (AB)  i3  defined  by  p  =1.  We  assume  that  the  duct  itself  is 
fitted  to  an  extremely  long  shock  tube  in  which,  by  means  of  an 
ideal  device,  the  flow  behind  the  shock  is  made  to  arrive  to  the 
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exit  without  a  sizeable  boundary  layer  and  practically  uniform 


across,  save  for  the  vanishing  of  the  velocity  on  the  wall.  At 
t=0,  the  shock  (which  we  will  call  the  precursor  shock)  is  just 
out  of  the  shock  tube  and  into  the  duct.  At  every  value  of  t,  we 
will  define  its  position  as 


c  =  p  (8  ,t)  (104) 

and  we  will  compute  the  flow  in  the  duct  (that  is,  between 
p=1  and  p=c  and  between  8=6  and  8=ir/2)  assuming  that  the  super¬ 
sonic  flow  at  the  entrance  o?  the  duct  remains  unchanged. 

The  computational  variables  X,  Y  and  T  are  defined  in 
terms  of  p,  8  and  t  as  follows: 

it  .  tanhta  ( Y— 1  )  ]  it 

8  =  (~  -  8  ) -  +  ~ 

2  o  tanh  a  2 

p  =  1  +  [c(Y, T)  -  1]X  (105) 

t  =  T 
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The  calculation  is  started  at  a  small,  positive  value  of 
t,  when  the  precursor  shock  has  already  moved  from  AB.  Initial¬ 
ly,  the  shock  is  assumed  to  lie  on  a  o=constant  line  (with  the 

value  of  the  constant,  c  ,  slightly  greater  than  1),  and  all 

o 

parameters  pertinent  to  the  shock  are  assumed  equal  to  their 
values  at  t=0  (since  the  flow  behind  the  shock  is  uniform,  the 
shock  Mach  number  defines  pressure,  velocity  and  entropy).  Uni¬ 
form  flow  is  assumed  between  AB  and  the  shock,  if  the  gas  is 
inviscid.  The  velocity  component,  v  is  made  equal  to  zero 
throughout.  Initially,  the  computational  region  is  divided  into 
two  strips  only  along  the  X-axis.  Along  the  Y-axis  we  consider 
as  many  partitions  as  necessary  to  provide  sufficient  resolution. 
Since  the  flow  is  viscous,  we  need  some  modifications  to  the  ini¬ 
tial  conditions  near  the  wall,  to  account  for  the  vanishing  of 
the  velocity  at  the  wall.  The  precursor  shock  cannot  reach  the 
wall;  the  perturbation  front,  elsewhere  in  the  form  of  a  shock, 
becomes  a  characteristic  at  the  wall,  somewhat  smeared  out  by 
viscous  diffusion.  Therefore,  on  the  initial  p=c  line  which 
represents  the  shock,  P  is  assumed  equal  to  zero  at  the  wall.  On 
the  next  wall  point,  P  is  taken  equal  to  one  half  of  its  value 
behind  the  shock.  The  wall  temperature  is  assumed  equal  to  1  at 
any  time.  The  entropy  is  defined  accordingly.  The  u-velocity 
component  is  taken  equal  to  one  half  of  its  value  behind  the 
shock  along  the  second  9=constant  line  from  the  wall.  The  calcu¬ 
lation  proceeds  as  detailed  above.  We  note  that  the  centerline 
is  computed  as  a  symmetry  line,  not  as  a  rigid  wall.  For  computa¬ 
tional  purposes,  the  3hock  geometry  is  prolonged  to  the  wall  with 
a  constant  value  of  c,  and  the  wall  point  is  computed  as  any  in¬ 
terior,  shockless  point,  taking  advantage  of  the  fact  that  the 
state  of  the  gas  in  front  of  the  perturbed  region  is  known.  The 
same  procedure  is  automatically  applied  to  any  other  point  on  the 
perturbation  front,  if  the  shock  happens  to  lose  its  strength 
completely. 

As  the  calculation  proceeds,  the  number  of  grid  intervals 
in  the  X-direction  is  doubled  every  time  (c-1)  on  the  centerline 
exceeds  1.4  times  its  initial  value  or  the  value  it  had  at  the 
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previous  doubling,  until  the  total  number  of  intervals  is  16. 

Certain  features  of  the  flow  are  common  to  the  inviscid 
and  viscous  models.  An  expansion  appears  from  the  beginning  near 
the  inlet  of  the  duct;  the  region  of  maximum  expansion  moves  from 
left  to  right,  but  at  a  slower  speed  than  the  precursor  shock. 
Consequently,  even  with  the  precursor  shock  losing  strength,  the 
pressure  behind  it  remains  higher  than  the  lowest  pressure  al¬ 
ready  attained  along  the  duct.  The  particles  are  accelerated  and 
then  decelerated  again.  A  compression  wave  appears,  which  tends 
to  steepen  up,  as  every  compression  wave  does,  and  another  shock 
results  eventually. 

In  the  presence  of  viscosity,  the  reccmpression  wave  and 
the  secondary  3hock  strongly  interact  with  the  boundary  layer  in 
the  process  of  formation.  The  latter  thickens  and  separates  very 
soon.  Between  the  main  stream  and  the  wall,  a  wide  dead-water 
region  appears,  where  the  pressure  tends  to  equalize  the  ambient 
pressure  (in  front  of  the  precursor  shock).  When  the  separated 
flow  becomes  steady,  it  is  what  is  commonly  defined  as  a  plume. 
From  the  separation  point  on,  the  plume  is  insensitive  to  the 
wall  geometry  and  the  flow  inside  it  is  essentially  inviscid. 
The  results  of  the  calculation  of  a  steady,  inviscid  flow  in  a 
plume  can  be  used  to  judge  whether  the  viscous  calculation  ap¬ 
proaches  its  theoretical  asymptote,  and  how  well. 

In  the  present  example,  a=2,  and  the  sector  between  9  =  9 
and  9=it/2  is  divided  into  16  intervals.  Consequently,  we  obtain 
a  fair  accumulation  of  grid  lines  near  the  wall,  and  still  work 
with  a  reasonably  small  number  of  lines.  It  i3  clear,  however, 
that  the  resolution  is  well  below  the  limits  which  are  usually 
recommended  for  a  good  description  of  Reynolds  number  effects; 
the  lack  of  an  adequate  resolution  is  particularly  felt  in  the 
vicinity  of  the  plume  shear  layer.  Nevertheless,  the  present 
results  are  very  encouraging,  just  because  very  good  qualitative 
results  are  obtained  with  such  a  coarse  mesh. 

Based  on  realistic  values  of  the  inner  diameter  of  the  shock  tube 

(25.4  mm),  of  the  kinematic  viscosity  of  air  (.0016  m  /sec)  and 

of  the  speed  of  sound  of  the  gas  at  rest  (360  m/sec) ,  R  turns 

e 

out  to  be  of  the  order  of  2500.  For  such  a  Reynolds  number,  the 
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resolution  of  the  grid  near  the  wall  cannot  provide  accurate  de¬ 
tails  of  the  boundary  layer,  such  as  needed,  for  example,  to  com¬ 
pute  the  skin  friction,  but  it  is  sufficient  to  furnish  an  ade¬ 
quate  picture  of  the  layer  in  the  general  context  of  the  flow. 

The  Prandtl  number  is  taken  equal  to  1,  and  y  equal  to 

1.4.  The  values  of  r  and  9  are  2.2603  and  1.2,  respectively. 

o  o 

The  temperature  at  the  wall  is  assumed  to  be  equal  to  1,  that  is, 
to  the  temperature  in  the  gas  at  rest.  The  value  of  P  at  the  in¬ 
let  is  1.6487.  Consistent  values  of  u  and  S  are  u=1.6551, 
S=0.1702.  The  Mach  number  at  the  inlet  is  barely  supersonic, 
M=1.040.  The  flow  is  assumed  to  be  two-dimensional. 

The  evolution  of  the  flow  i's  shown  in  Figs.  2  through  13. 

In  Figs.  2  through  8,  isobars  are  plotted.  The  left  vertical 

line  is  the  inlet  cross-section,  the  lower  curved  line  is  the 
wall  of  the  duct,  the  upper  horizontal  line  is  the  duct  center- 
line  (on  which  notches  indicate  multiples  of  the  unit  length); 
and  the  right  boundary  is  the  precursor  shock.  The  line  marked 
with  0  is  the  sonic  line.  Step  800  is  typical  of  the  first  phase 
of  evolution,  showing  an  accented  minimum  of  pressure  on  the 
wall;  velocity  vectors  drawn  for  that  step  indicate  the  beginning 

of  recirculation  in  the  boundary  layer  around  the  region  of 

minimum  pressure  (Fig.  9).  At  step  1200,  the  minimum  pressure 
has  moved  to  the  centerline,  and  a  steepening  up  of  the  pressure, 
all  across  the  duct,  is  evident.  At  step  1400,  an  imbedded  shock 
is  fitted  (marked  in  the  figure  by  +  ).  The  initial  fitting  of 
the  shock  is  rather  arbitrary,  but  such  arbitrariness  is  not  res¬ 
trictive.  In  general,  an  imbedded  shock  is  fitted  on  any 
9=constant  line,  in  the  middle  of  an  interval  where  the  differ¬ 
ence  in  P  exceeds  0.6,  as  it  is  suppressed  if  its  normal  Mach 
number  becomes  less  than  1  or  if  locally  the  shock  stretches  over 
more  than  two  X-intervals  over  a  single  Y-interval. 

Note  that  the  computational  technique  allows  the  shock  to 
end  inside  the  flow  field,  and  that  isobars  between  the  end-point 
of  the  shock  and  the  wall  shape  up  in  a  way  which  is,  at  least 
qualitatively,  typical  of  shock-boundary-layer  interactions. 
Velocity  vectors  at  this  3tep  (Fig.  10)  show  a  pronounced  recir- 
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culation  bubble. 

At  step  2400,  the  flow  field  is  mostly  subsonic;  the 
supersonic  flow  is  confined  to  a  jet-like  region  (see  Fig.  11, 
where  lines  of  constant  Mach  number  are  shown,  and  the  veocity 
vectors  in  Fig.  12). 

At  step  3400  the  loss  of  accuracy  consequent  to  loss  of 
resolution  is  evident;  the  jet  region,  where  the  imbedded  shock 
is  still  present,  is  covered  by  three  mesh  intervals  only.  The 
Mach  number  distribution  (Fig.  13)  has  the  correct  trend,  but  the 
transversal  gradient  is  spread  out  too  widely.  Nevertheless, 
note  the  formation  of  a  transversal  pressure  gradient,  very  simi¬ 
lar  to  the  pattern  found  in  the  isobars  computed  for  a  steady 
plume  originated  by  the  same  duct,  with  separation  taking  place 
exactly  where  the  pressure  at  the  wall  equals  the  ambient  pres¬ 
sure  (Fig.  14). 

To  conclude,  we  show  an  impressive  picture  of  streamlines 
at  this  step  in  Fig.  15. 
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